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Abstract 

Calibration of radio interferometric observations becomes increasingly difficult towards lower frequencies. Below ~ 300 
MHz, spatially variant refractions and propagation delays of radio waves traveling through the ionosphere cause phase 
, rotations that can vary significantly with time, viewing direction and antenna location. In this article we present a 

^ I' description and first results of SPAM (Source Peeling and Atmospheric Modeling), a new calibration method that 

attempts to iteratively solve and correct for ionospheric phase errors. To model the ionosphere, we construct a time- 
variant, 2-dimensional phase screen at fixed height above the Earth's surface. Spatial variations are described by a 
truncated set of discrete Karhunen-Loeve base functions, optimized for an assumed power-law spectral density of free 
^ . electrons density fiuctuations, and a given configuration of calibrator sources and antenna locations. The model is 

constrained using antenna-based gain phases from individual self-calibrations on the available bright sources in the 
field-of-view. Application of SPAM on three test cases, a simulated visibility data set and two selected 74 MHz VLA 
, data sets, yields significant improvements in image background noise (5-75 percent reduction) and source peak fluxes 

^ ■ (up to 25 percent increase) as compared to the existing self-calibration and field-based calibration methods, which 

1^-^ ' indicates a significant improvement in ionospheric phase calibration accuracy. 
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">) ' 1. Introduction Ionospheric effects on LF interferometric observations 

. have usually been ignored for several reasons: (i) the res- 

a> . Radio waves of cosmic origin are influenced by the Earth's ^^^^-^^^ sensitivity of the existing arrays were gener- 

; atmosphere before detection at ground level. At low fre- ^jjy to be affected, (ii) existing calibration algo- 

- quencies (LF; < 300 MHz), the dommant effects are refrac- ^^^-^^^ self-calibration) appeared to give reasonable 

tion, propagation delay and Faraday rotation caused^ ^gs^l^s ^^os^ of ^j^g time, and (iii) a lack of computing power 

the ionosphere (e.g., Thompson, Moran & Swenson [20011 ^^^^ ^he needed calculations prohibitly expensive. During 

TMS2001 from here on). For a ground-based mterferom- ^j^g j^g^ 15 years, two large and more sensitive LF arrays 

eter (array from here on) observing a LF cosmic source, ^^^ve become operational: the VLA at 74 MHz (Kassim et 

the ionosphere is the mam source of phase errors m the [Wf|) and the GMRT at 153 and 235 MHz (Swarup 

visibilities. Amplitude errors may also arise under severe Observations with these arrays have demonstrated 

ionospheric conditionsdue to diffraction or focussing (e.g., ^hat ionospheric phase errors arc one of the main limiting 

Jacobson & Erickson [1992]) . ^ ^ factors for reaching the theoretical image noise level. 

The ionosphere causes propagation delay differences be- „ ,.i r r,i ir, i 

, , 1, • • 1 • ii • • ror optimal performance oi these and future large ar- 

tween array elements, resulting m phase errors m the visi- \ „ ^, , , „„. „ t^htx ^ oV.^ x \ 

u-i-i- rn, J , T 2 / ^ r 1 \ rays with LF capabihties (such as LOFAR, LWA and SKA), 

bilities. 1 he delay per array element (antenna from here on) . • i-i :• i -^i l^ ^ / 

J J ii 1- r • 1 J. /T o\ ii 1 ii • 1 it IS crucial to use calibration algorithms that can properly 

depends on the hne-oi-sight (LoS) through the ionosphere, , , , . i ■ ^ •, r ^^ ■ ■ 

\ ^, J. , . ' , . . 1- model and remove ionospheric contributions irom the visi- 

and thereiore on antenna position and viewing direction. , i i i i i-r ^- /r-i ^ ^ i aaa j k • ,i 

m, i-i f? u ■ 1 bihties. lield-based calibration (Cotton et al. I2004D is the 

Ihe calibration oi LI observations requires phase correc- . , . ,. . , . i-, ;. » • ■ .^ ^ 

.. , j-u c ij -(• • ^t?\t\ r 1, i single existing ionospheric calibration & imaging method 

tions that vary over the neld-oi- view (i'oV) ot each antenna. ^, ^ . °^ ^■ ^- ^ ^ ^ ^ ^■^ ^■ 

^ vi i- ii 1 ii i 1 i • • i 1 that incorporates direction-dependent phase calibration. 

Calibration methods that determine lust one phase correc- • ^ , . , , r n i • i ^ -trr a 

r r „ r , , n-i ir 1-1 i- 1 his tcchmquc has been succcstully applied to many VLA 

tion tor the mil 1*0 V of each antenna (like selt-calibration; ~a ^^tt ^^ ^ t ^ ■ ^■ ^ ^ ^ ■ r -^i 

r) p r> JT, J H A6^ h i.u c • fc 74 MHz data sets, but is limited by design lor use with 

e.g., see Pearson & Readhead 119841) are thereiore msulh- , i J j t> 

. ' ' relatively compact arrays. 

cient. J f J 

In Section [21 we discuss ionospheric calibration in more 

Send offprint requests to: H.T. Intema, detail. In Section [3l we present a detailed description of 

e-mail: intema@strw.leidenuniv.nl SPAM, a new ionospheric calibration method that is appli- 
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cable to LF observations with relatively larger arrays. In 
Section 21 we present the first results of SPAM calibration 
on simulated and real VLA 74 MHz observations and com- 
pare these with results from self-calibration and field-based 
calibration. A discussion and conclusions are presented in 
Section [5l 

2. Ionosphere and Calibration 

In this Section, we describe some physical properties of the 
ionosphere, the phase effects on radio interferometric obser- 
vations and requirements for ionospheric phase calibration. 



2.1. The Ionosphere 

The ionosphere is a partially ionised layer of gas between 
~ 50 and 1000 km altitude over the Earth's surface (e.g., 
Davies I1990p . It is a dynamic, inhomogeneous medium, 
with electron density varying as a function of position and 
time. The state of ionization is mainly influenced by the 
Sun through photo-ionization at UV and short X-ray wave- 
lengths and through injection of charged particles from the 
solar wind. Ionization during the day is balanced by recom- 
bination at night. The peak of the free electron density is 
located at a height around 300 km. The free electron column 
density along a LoS through the ionosphere is generally re- 
ferred to as total electron content^ or TEC. The TEC unit 
(TECU) is 10^^ m~^ which is a typically observed value at 
zenith during nighttime. 

The refraction and propagation delay are caused by a 
varying refractive index n of the ionospheric plasma along 
the wave trajectory. For a cold, coUisionless plasma without 
magnetic field, n is a function of the free electron density 
rig and is defined by (e.g., TMS2001) 



(1) 



with V the radio frequency and the plasma frequency. 



given by 



e_ I 

2tt Y £9"^' 



(2) 



with e the electron charge, m the electron mass, Cq the vac- 
uum permittivity. Typically, for the ionosphere, ranges 
from 1-10 MHz, but may locally rise up to ~ 200 MHz in 
the presence of sporadic E-laycrs (clouds of unusually high 
free electron density) . Cosmic radio waves with frequencies 
below the plasma frequency are reflected by the ionosphere 
and do not reach the Earth's surface. For higher frequencies, 
the spatial variations in electron density cause local refrac- 
tions of the wave (Snell's Law) as it travels through the 
ionosphere, thereby modifying the wave's trajectory. The 
total propagation delay, integrated along the LoS, results 
in a phase rotation given by 



-^/(n-l)d^, 



(3) 



with c the speed of light in vacuum. For frequencies v ^ v^^ 
this can be approximated by 



where the integral over on the right is the TEC along 
the LoS. Note that this integral depends on the wave's tra- 
jectory, and therefore on local refraction. Because the re- 
fractive index is frequency-dependent, the wave's trajectory 
changes with frequency. As a consequence, the apparent 
scaling relation 0'°'^ oc i^^^ from Equation[4]is only valid to 
first order in frequency. 

Although bulk changes in the large scale TEC (e.g., 
a factor of 10 increase during sunrise) have the largest 
amplitudes, the fluctuations on relatively small spatial 
scales and short temporal scales are most troublesome for 
LF interferometric observations. Most prominent are the 
traveling ionospheric disturbances (TIDs), a response to 
acoustic-gravity waves in the neutral atmosphere (e.g., van 
Velthoven I1990p . Typically, medium-scale TIDs are ob- 
served at heights between 200 and 400 km, have wave- 
lengths between 250 and 400 km, travel with near- 
horizontal velocities between 300 and 700 km h~^ in 
any direction and cause 1-5 percent variations in TEC 
(TMS2001). 

The physics behind fluctuations on the shortest spa- 
tial and temporal scales is less well understood. Temporal 
and spatial behaviour may be coupled through quasi-frozen 
patterns that move over the area of interest with a certain 
velocity and direction (Jacobson & Erickson ll992p . Typical 
variations in TEC are on the order of 0.1 percent, observed 
on spatial scales of tens of kilometers down to a few km, 
and time scales of minutes down to a few tens of seconds. 
The statistical behaviour of radio waves passing through 
this medium suggests the presence of a turbulent layer with 
a power-law spectral density of free electron density fluc- 
tuations P„^((?) oc q-"' (e.g., TMS2001), with q = \q\ the 
magnitude of the 3-dimensional spatial frequency. P„ (g) is 
defined in units of electron density squared per spatial fre- 
quency. The related 2-dimensional structure function of the 
phase rotation 4> of emerging radio waves from a turbulent 
ionospheric layer is given by 



(5) 



— / v'dl = 



AL 



(4) 



where x and x -\- r are Earth positions, r = \r\ is the hori- 
zontal distance between these two points, (. . .) denotes the 
expected value and 7 = a — 2. For pure Kolmogorov tur- 
bulence, a = 11/3, therefore 7 = 5/3. 

Using differential Doppler-shift measurements of satel- 
lite signals, van Velthoven p990p found a power-law rela- 
tion between spectral amplitude of small-scale ionospheric 
fluctuations and latitudinal wave-number with exponent 
a/2 — 3/2. Combining with radio interferometric obser- 
vations of apparent cosmic source shifts, van Velthoven 
derived a mean height for the ionospheric perturbations 
of 200-250 km. Through analysis of differential apparent 
movement of pairs of cosmic sources in the VLSS, Cohen & 
Rottgering (|2008p find typical values for 7/2 of 0.50 during 
nighttime and 0.69 during daytime. Direct measurement of 
phase structure functions from different GPS satellites (van 
der Tol, unpublished) shows a wide distribution of values for 
7 that peaks at ~ 1.5. On average, these results indicate 
the presence of a turbulent layer below the peak in the free 
electron density that has more power in the smaller scale 
fluctuations than in the case of pure Kolmogorov turbu- 
lence. Note that for individual observing times and loca- 
tions, the behaviour of small-scale ionospheric fluctuations 
may differ significantly from this average. 
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2.2. Image Plane Effects 

Interferometry uses the phase differences as measured on 
basehnes to determine the angle of incident waves, and is 
therefore only sensitive to TEC differences. A baseline is 
sensitive to TEC fluctuations with linear sizes that arc com- 
parable to or smaller than the baseline length. At 75 MHz, a 
0.01 TECU difference on a baseline causes a ~ I radian visi- 
bility phase error (Equation 2]). Because the observed TEC 
varies with time, antenna position and viewing direction, 
visibility phases are distorted by time-varying differential 
ionospheric phase rotations. 

An instantaneous spatial phase gradient over the array 
in the direction of a source causes an apparent position shift 
in the image plane (e.g., Cohen & R5ttgering |2008p . but no 
source deformation. If the spatial phase behaviour deviates 
from a gradient, this will also distort the apparent shape of 
the source. Combining visibilities with different time labels 
while imaging causes the image plane effects to be time- 
averaged. A non-zero time average of the phase gradient 
results in a source shift in the final image. Both a zero- 
mean time variable phase gradient and higher order phase 
effects cause smearing and deformation of the source image, 
and consequently a reduction of the source peak fiux (see 
Cotton & Condon [20021 for an example). In the latter case, 
if the combined phase errors behave like Gaussian random 
variables, a point source in the resulting image experiences 
an increase of the source width and reduction of the source 
peak flux, but the total flux (the integral under the source 
shape) is conserved. 

For unresolved sources, the Strehl ratio is defined as the 
ratio of observed peak flux over true peak flux. In case of 
Gaussian random phase errors, the Strehl ratio R is related 
to the RMS phase error by (Cotton et al. [MM)) 

i? = exp(^-^^. (6) 

A larger peak flux is equivalent to a smaller RMS phase 
error. This statement is more generally true, because all 
phase errors cause scattering of source power into sidelobes. 

A change in the apparent source shape due to iono- 
spheric phase errors leads to an increase in residual side- 
lobes after deconvolution. Deconvolution subtracts a time- 
averaged source image model from the visibility data at 
all time stamps. In the presence of time- variable phase er- 
rors, the mean source model deviates from the apparent, 
instantaneous sky emission and subtraction is incomplete. 
Residual sidelobes increase the RMS background noise level 
and, due to its non-Gaussian character, introduce struc- 
ture into the image that mimics real sky emission. In LF 
observations, due to the scaling relation of the dirty beam 
with frequency (width oc v~^), residual sidelobes around 
bright sources can be visible at significant distances from 
the source. 

2.3. Ionospheric Phase Calibration 

Lonsdale pOOSp discussed four different regimes for (instan- 
taneous) ionospheric phase calibration, depending on the 
different linear spatial scales involved. These scales are the 
array size the scale size S of ionospheric phase fluctu- 
ations and the projected size V of the field-of-view (FoV) 
at a typical ionospheric height. We use the term compact 



array when A <^ S and extended array when A > 5. Note 
that these definitions change with ionospheric conditions, 
so there is no fixed linear scale that defines the difference 
between compact and extended. A schematic overview of 
the different regimes is given in Figure |TJ 

The combination AV/ S"^ is a measure of the complexity 
of ionospheric phase calibration. Both S and V depend on 
the observing frequency v. For a power-law spectral den- 
sity of free electron density fluctuations (see Section 12. ip 
S scales with v , and for a fixed circular antenna aperture 
V scales with . Therefore, AV/ S'^ scales with , sig- 
nalling a rapid increase in calibration problems towards low 
frequencies. 

Under isoplanatic conditions (V^ <C 5), the ionospheric 
phase error per antenna does not vary with viewing direc- 
tion within the FoV, for both compact and large arrays 
(Lonsdale regimes 1 and 2, respectively). Phase-only self- 
calibration on short enough time-scales is sufficient to re- 
move the ionospheric phase errors from the visibilities. 

Under anisoplanatic conditions {V > S*), the iono- 
spheric phase error varies over the FoV of each antenna. A 
single phase correction per antenna is no longer sufficient. 
Self-calibration may still converge, but the resulting phase 
correction per antenna is a flux-weighted average of iono- 
spheric phases across the FoV (see Section [3T|) . Accurate 
self-calibration and imaging of individual very bright and 
relatively compact sources is therefore possible, even with 
extended arrays (see Gizani, Cohen & Kassim [20051 for an 
example). For a compact array (Lonsdale regime 3), the 
FoV of different antennas effectively overlap at ionospheric 
height. The LoS of different antennas towards one source 
run close and parallel through the ionosphere. For an ex- 
tended array (Lonsdale regime 4), the FoV of different an- 
tennas may partially overlap at ionospheric height, but not 
necessarily. Individual LoS from widespread antennas to 
one source may trace very different paths through the iono- 
sphere 

In regime 3, ionospheric phases behave as a spatial gra- 
dient over the array that varies with viewing direction. 
This causes the apparent position of sources to change 
with time and viewing direction, but no source deforma- 
tion takes place. The 3-dimensional phase structure of the 
ionosphere can be effectively reduced to a 2-dimensional 
phase screen, by integrating the free electron density along 
the LoS (Equation |4|) . Radio waves that pass the virtual 
screen experience an instantaneous ionospheric phase rota- 
tion depending on the pierce point position (where the LoS 
pierces the phase screen). When assuming a fixed number 
of required ionospheric parameters per unit area of phase 
screen, calibration of a compact array requires a minimal 
number of parameters because each antenna illuminates the 
same part of the phase screen. 

In regime 4, the dependence of ionospheric phase on 
antenna position and viewing direction is more complex. 
This causes source position shifts and source shape defor- 
mations that both vary with time and viewing direction. 
A 2-dimcnsional phase screen model may still be used, but 
only when the dominant phase ffuctuations originate from 
a restricted height range Ah <C 5 in the ionosphere. The 
concept of a thin layer at a given height is attractive, be- 
cause it reduces the complexity of the calibration problem 
drastically. When using an airmass function to incorporate 
a zenith angle dependence, the spatial phase function is in 
effect reduced to 2 spatial dimensions. Generally, a phase 
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Figurcl: Schematic overview of the difFcrcnt cahbration regimes as discussed by Lonsdale (|2005p . For clarity, only two 
spatial dimensions and one calibration time interval are considered. In this overview, the array is represented by three 
antennas at ground level, looking through the ionospheric electron density structure (grey bubbles) with individual ficlds- 
of-view (red, green and blue areas). Due to the relatively narrow primary beam patterns in regimes 1 and 2 (top left and 
top right, respectively), each individual antenna 'sees' an approximately constant TEC across the FoV. The relatively wide 
primary beam patterns in regimes 3 and 4 [bottom left and bottom right, respectively) causes the antennas to 'see' TEC 
variations across the FoV. For the relatively compact array configurations in regimes 1 and 3, the TEC variation across 
the array for a single viewing direction within the FoV is approximately a gradient. For the relatively extended array 
configurations in regimes 2 and 4, the TEC variation across the array for a single viewing direction differs significantly 
from a gradient. The consequences for calibration of the array are discussed in the text. 



screen in regime 4 requires a larger number of model pa- 
rameters than in regime 3, because the phase screen area 
illuminated by the total array is larger. 

It is currently unclear under which conditions a 2- 
dimensional phase screen model becomes too inaccurate to 
model the ionosphere in regime 4. For very long baselines 
or very severe ionospheric conditions, a full 3-dimensional 
ionospheric phase model may be required, where iono- 
spheric phase corrections need to determined by ray- 
tracing. Such a model is likely to require many more pa- 
rameters than can be extracted from radio observations 
alone. To first order, it may be sufhcient to extend the 
phase screen model with some form of height-dependence. 
Examples of such extensions are the use of several phase 
screens at different heights (Anderson I2006P or introducing 
smoothly varying partial derivatives of TEC or phase as a 
function of zenith angle (Noordam et al. \in preparation^ . 

Calibration needs to determine corrections on suffi- 
ciently short time scales to track the ionospheric phase 



changes. The phase rate of change depends on the intrin- 
sic time variability of the TEC along a given LoS and on 
the speed of the LoS from the array antennas through the 
ionosphere while tracking a cosmic source. The latter may 
range up to ~ 100 km h^^ at 200 km height. The exact 
requirements on the time resolution of the calibration are 
yet to be determined. In principle, the time-variable iono- 
spheric phase distortions needs to be sampled at least at 
the Nyquist frequency. However, during phase variations 
of large amplitude (^ 1 radian), 2tt radian phase winding 
introduces periodicity on much shorter time scales. To suc- 
ccsfully unwrap phase winds, at least two corrections per 
27r radian phase change are required. 

2.4. Proposed and Existing Ionospheric Calibration Schemes 

Schwab (|1984p and Subrahmanya (|199ip have proposed 
modifications to the self-calibration algorithm to support 
direction-dependent phase calibration. Both methods dis- 
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cuss the use of a spatial grid of interpolation nodes (addi- 
tional free parameters) to characterize the spatial variabil- 
ity of the ionospheric phase error. Schwab suggests to use 
a different set of nodes per antenna, while Subrahmanya 
suggests to combine these sets by positioning them in a 
quasi-physical layer at fixed height above the Earth's sur- 
face (this to reduce the number of required nodes when the 
FoVs from different antennas overlap at ionospheric height). 
Neither of both proposed methods have been implemented. 

Designed to operate in Lonsdale regime 3, field-based 
calibration by Cotton et al. (2004j) is the single existing 
implementation of a direction-dependent ionospheric phase 
calibration algorithm. Typically, for each time interval of 
1-2 minutes of VLA 74 MHz data, the method measures 
and converts the apparent position shift of 5-10 detectable 
bright sources within the FoV into ionospheric phase gradi- 
ents over the array. To predict phase gradients in arbitrary 
viewing directions for imaging of the full FoV, an indepen- 
dent phase screen per time interval is fitted to the measured 
phase gradients. The phase screen is described by a 5 term 
basis of Zernike polynomials (up to second order, excluding 
the constant zero order). 

Field-based calibration has been used to calibrate 
74 MHz VLA observations, mostly in B-configuration (e.g., 
Cohen et al. 12007) but also several in A-configuration (e.g., 
Cohen et al. i2003,, 2004]) . Image plane comparison of field- 
based calibration against self-calibration shows an overall 
increase of source peak fiuxes (in some cases up to a factor 
of two) and reduction of residual sidelobes around bright 
sources, a clear indication of improved phase calibration 
over the FoV (Cotton & Condon [2002p . The improved over- 
all calibration performance sometimes compromises the cal- 
ibration towards the brightest source. 

Zernike polynomials are often used to describe abbera- 
tions in optical systems, because lower order terms match 
well with several different types of wavefront distortions, 
and the functions are an orthogonal set on the circular do- 
main of the telescope pupil. Using Zernike polynomials to 
describe an ionospheric phase screen may be less suitable, 
because they are not orthogonal on the discrete domain of 
pierce points, diverge when moving away from the field cen- 
ter and have no relation to ionospheric image abberations 
(except for first order, which can model a large scale TEC 
gradient). Non-orthogonality leads to interdependence be- 
tween model parameters, while divergence is clearly non- 
physical and leads to undesirable extrapolation properties. 

For extended LF arrays or more severe ionospheric con- 
ditions, the ionospheric phase behaviour over the array for 
a given viewing direction is no longer a simple gradient. 
Under these conditions, performance of field-based calibra- 
tion degrades. For the 74 MHz VLA Low-frequency Sky 
Survey (VLSS; Cohen et al. I2007p . field-based calibration 
was unable to calibrate the VLA in B-configuration for 
about 10-20% of the observing time due to severe iono- 
spheric conditions. Observing at 74 MHz with the ~ 3 times 
larger VLA A-configuration leads to a relative increase in 
the failure rate of field-based calibration. This is to be ex- 
pected, as the larger array size results in an increased prob- 
ability for the observations to reside in Lonsdale regime 4. 

The presence of higher order phase structure over the 
array in the direction of a calibrator requires an antenna- 
based phase calibration rather than a source position shift 
to measure ionospheric phases. The calibration methods 
proposed by Schwab and Subrahmanya (see above) do al- 



low for higher order phase corrections over the array and 
could, in principle, handle more severe ionospheric condi- 
tions. An alternative approach is to use the peeling tech- 
nique (Noordam I2004p . which consist of sequential self- 
calibrations on individual bright sources in the FoV. This 
yields per source a set of time- variable antenna-based phase 
corrections and a source model. Because the peeling correc- 
tions are applicable to a limited set of viewing directions, 
they need to be interpolated in some intelligent way to arbi- 
trary viewing directions while imaging the full FoV. Peeling 
is described in more detail in Section [3.31 

Noordam (|2004p has proposed a 'generalized' self- 
calibration method for LOFAR (e.g., Rottgering et al. l2006p 
that includes calibration of higher order ionospheric phase 
distortions. Similar to 'classical' self-calibration, instru- 
mental and environmental (including ionospheric) parame- 
ters are estimated by calibration against a sky brightness 
model. Sky model and calibration parameters are itera- 
tively updated to converge to some final result. Uniqueness 
of the calibration solution is controlled by putting re- 
strictions on the time-, space- and frequency behaviour 
of the fitted parameters. The effects of the ionosphere 
are modeled in a Minimum Ionospheric Model (MIM; 
Noordam |m preparation^ , which is yet to be defined in de- 
tail. The philosophy of the MIM is to use a minimal number 
of physical assumptions and free parameters to accurately 
reproduce the observed effects of the ionosphere on the vis- 
ibilities for a wide-as-possible range of ionospheric condi- 
tions. The initial MIM is to be constrained using peeling 
corrections. 

3. Method 

SPAM, an abbreviation of 'Source Peeling and Atmospheric 
Modeling', is the implementation of a new ionospheric cali- 
bration method, combining several concepts from proposed 
and existing calibration methods. SPAM is designed to op- 
erate in Lonsdale regime 4 and can therefore also operate 
in regimes 1 to 3. It uses the calibration phases from peel- 
ing sources in the FoV to constrain an ionospheric phase 
screen model. The phase screen mimics a thin turbulent 
layer at a fixed height above the Earth's surface, in con- 
cordance with the observations of ionospheric small-scale 
structure fScction l2.ip . The main motivation for this work 
was to test several aspects of ionospheric calibration on 
existing VLA and GMRT data sets on viability and quali- 
tative performance, and thereby support the development 
of more advanced calibration algorithms for future instru- 
ments such as LOFAR. 

Generally, the instantaneous ionosphere can only be 
sparsely sampled, due to the non-uniform sky distribution 
of a limited number of suitable calibrators and an array 
layout that is optimized for UV-coverage rather than iono- 
spheric calibration. To minimize the error while interpolat- 
ing to unsampled regions, an optimal choice of base func- 
tions for the description of the phase screen is of great im- 
portance. Based on the work by van der Tol & van der Veen 
(|2007p . we use the discrete Karhunen-Loeve (KL) transform 
to determine an optimal set of base 'functions' to describe 
our phase screen. For a given pierce point layout and an 
assumed power-law slope for the spatial structure function 
of ionospheric phase fluctuations (see Section 12. ip , the KL 
transform yields a set of base vectors with several impor- 
tant properties: (i) the vectors are orthogonal on the pierce 
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point domain, (ii) truncation of the set (reduction of the 
model order) gives a minimal loss of information, (iii) in- 
terpolation to arbitrary pierce point locations obeys the 
phase structure function, and (iv) spatial phase variabil- 
ity scales with pierce point density, i.e., most phase screen 
structure is present in the vicinity of pierce points, while it 
converges to zero at infinite distance. More detail on this 
phase screen model is given in Section [3.41 

Because the required calibration time resolution is still 
an open issue, and the SPAM model does not incorporate 
any restrictions on temporal behaviour, independent phase 
screens are determined at the highest possible time resolu- 
tion (which is the visibility integration time resolution). 

SPAM calibration can be separated in a number of func- 
tional steps, each of which is discussed in detail in the sec- 
tions to follow. The required input is a spectral-mode visi- 
bility data set that has flux calibration and bandpass cali- 
bration applied, and radio frequency interference (RFI) ex- 
cised (see Lazio, Kassim & Perlev l^OOSI or Cohen et al. 120071 
for details). The SPAM recipe consists of the following 
steps: 

1. Obtain and apply instrumental calibration corrections 
for phase (Section l3.ip . 

2. Obtain an initial model of the apparent sky, to- 
gether with an initial ionospheric phase calibration 
(Section [321) • 

3. Subtract the sky model from the visibility data while 
applying the phase calibration. Peel apparently bright 
sources (Section 13. 3|) . 

4. Fit an ionospheric phase screen model to the peeling 
solutions (Section [3l4|) . 

5. Apply the model phases on a facet-to- facet basis during 
re-imaging of the apparent sky (Section [23]). 

Steps[3|to[5|define the SPAM calibration cycle, as the image 
produced in step [5] can serve as an improved model of the 
apparent sky in step[3l 

The scope of applications for SPAM is limited by a num- 
ber of assumptions that were made to simplify the current 
implementation: 

— The ionospheric inhomogeneities that cause significant 
phase distortions are located in a single, relatively nar- 
row height range. 

— There exists a finitely small angular patch size, which 
can be much smaller than the FoV of an individual 
antenna, over which the ionospheric phase contribu- 
tion is effectively constant. Moving from one patch to 
neighbouring patches results in small phase transitions 
(^ 1 radian). 

— There exists a finitely small time range, larger than the 
integration time interval of an observation, over which 
the apparent ionospheric phase change for any of the 
array antennas along any line-of-sight is much smaller 
than a radian. 

— The bandwidth of the observations is small enough to 
be effectively monochromatic, so that the ionospheric 
dispersion of waves within the frequency band is negli- 
gible. 

— Within the given limitations on bandwidth and inte- 
gration time, the array is sensitive enough to detect at 
least a few (> 5) sources within the target FoV that 
may serve as phase calibrators. 



— The ionospheric conditions during the observing run 
are such that self-calibration is able to produce a good 
enough initial calibration and sky model to allow for 
peeling of multiple sources. This might not work under 
very bad ionospheric conditions, but for the applications 
presented in this article it proved to be sufficient. 

— After each calibration cycle (steps [3| to [5|) , the calibra- 
tion and sky model are equally or more accurate than 
the previous. This implies convergence to a best achiev- 
able image. 

— The instrumental amplitude and phase contributions to 
the visibilities, including the antenna power patterns 
projected onto the sky towards the target source, are 
constant over the duration of the observing run. 

SPAM does not attempt to model the effects of ionospheric 
Faraday rotation on polarization products, and is therefore 
only applicable to intensity measurements (stokes I). 

In our implementation we have focussed on function- 
ality rather than processing speed. In its current form, 
SPAM is capable of processing quite large offline data sets, 
but is not suitable for real-time processing as is required 
for LOFAR calibration. SPAM relies heavily on function- 
ality available in NRAO's Astronomical Image Processing 
System (AIPS; e.g.. Bridle & Greisen [TgM)) . It consists of a 
collection of Python scripts that accesses AIPS tasks, files 
and tables using the ParselTongue interface (Kettenis et 
al. I2006p . Two main reasons to use AIPS are its familiar- 
ity and proven robustness while serving a large group of 
users over a 30 year lifetime, and the quite natural way 
by which the ionospheric calibration method is combined 
with polyhedron imaging (Per ley [19893 Cornwell & Per ley 
I1992|) . SPAM uses a number of 3'''* party Python libraries, 
like scipy, numpy and matplotlib for math and matrix op- 
erations and plotting. For non-linear least squares fitting of 
ionospheric phase models, we have adopted a Levenberg- 
Marquardt solver (LM; e.g.. Press et al. I1992p based on 
IDE's MPFIT package (Markwardt [20081) . 

3.1. Instrumental Phase Calibration 

Each antenna in the array adds an instrumental phase offset 
to the recorded signal before correlation. At low frequencies, 
changes in the instrumental signal path length (e.g., due 
to temperature induced cable length differences) are very 
small compared to the wavelength, therefore instrumental 
phase offsets are generally stable over long time periods 
(hours to days). SPAM requires removal of the instrumen- 
tal phase offsets from the visibilities prior to ionospheric 
calibration. 

Instead of directly measuring the sky intensity /(^, m) 
as a function of viewing direction cosines (/, m), an interfer- 
ometer measures an approximate Fourier transform of the 
sky intensity. For a baseline consisting of antennas i and j, 
the perfect response to all visible sky emission for a single 
time instance and frequency is given by the measurement 
equation (ME) for visibilities (e.g., TMS2001): 

T/y=y y/(Z,m)e-2'^-'["-'+"-"+"'-("-i)]^^, (7) 

where J indicates the imaginary part of a complex number, 
n = ^y^ — P — m?^ and v^j are baseline coordinates in 
the UV plane (expressed in wavelengths) parallel to / and 
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m, respectively, and w.^j is the perpendicular baseline co- 
ordinate along the LoS towards the chosen celestial phase 
tracking center at {l,m) = (0,0). In practise, these mea- 
surements are modified with predominantly antenna-based 
complex gain factors that may vary with time, frequency, 
antenna position and viewing direction. This modifies the 
ME into 



a^{l, m) rn) 



Til ™\ -27rjl'u .;+u .m-l-u) .{n-l)l dm 

[{i,m)e L ij ' »j 'I , 



(8) 



Determination of the gain factors is generally referred to 
as calibration. When known, only gain factors that do not 
depend on viewing direction can be removed from the vis- 
ibility data prior to image reconstruction by applying the 
calibration: 



(9) 



This operation is generally not possible for gain factors that 
do depend on viewing direction, because these gain factors 
cannot be moved in front of the integral in Equation [51 
One may still choose to apply gain corrections for a single 
viewing direction (e.g. to image a particular source), but 
the accuracy of imaging and deconvolution of other visible 
sources will degrade when moving away from the selected 
viewing direction. A solution for wide-field imaging and de- 
convolving in the presence of direction-dependent gain fac- 
tors is discussed in Section 

The standard approach for instrumental phase cali- 
bration at higher frequencies is to repeatedly observe a 
bright (mostly unresolved) source during an observing run. 
Antenna-based gain phase corrections ~ a^^ are esti- 
mated by minimizing the weighted difference sum S be- 
tween observed visibilities Va^ and source model visibilities 



T^modcl ^ 

CALIB): 



(e.g., TMS2001; implemented in AIPS task 



i j>i 



del 



(10) 



with W^j the visibility weight (reciproke of the uncer- 
tainty in the visibility measurement), g^ = e"^* and p the 
power of the norm (typically 1 or 2). The source model 
visibilities V^™°^'^^ are calculated using Equation [7] with 

rn) = /™°'^°'(^, m). The phase corrections consist of 
an instrumental and an atmospheric part. The corrections 
are interpolated in time and applied to the target field vis- 
ibilities, under the assumptions that the instrumental and 
atmospheric phase offsets vary slowly in time, and that the 
atmospheric phase offsets in the direction of the target are 
equal to those in the direction of the calibrator. 

At low frequencies, there are two complicating factors 
for the standard approach: (i) the FoV around the calibra- 
tor source is large and includes many other sources, and 
(ii) the ionospheric phase offset per antenna changes signif- 
icantly with time and viewing direction. The former can be 
overcome by choosing a very bright calibrator source with 
a flux that dominates over the combined flux of all other 
visible sources on all baselines. For the VLSS (Cohen et al. 
I2007p . the 17,000 Jy of Cygnus A was more than sufficient 
to dominate over the total apparent flux of 400 — 500 Jy in a 



typical VLSS field. The latter requires filtering of the phase 
corrections to extract only the instrumental part, which is 
then applied to the target field visibilities. 

For SPAM, we have adopted an instrumental phase cal- 
ibration method that is very similar to the procedure used 
for field-based calibration (Cotton et al. I2004p . Antenna- 
based phase corrections are obtained on the highest possi- 
ble time resolution by calibration on a very bright source 
k using the robust LI norm (Equation [10] with p ~ 1; 
Schwab I198ip . A phase correction for antenna i at 
time interval n consist of several contributions: 



6^ = -^r'-^ 



lion 



rkn 



lambig 



(11) 



where the instrumental and ionospheric phase corrections. 



and 



respectively, are assumed to be constant 



resp. vary with time and antenna position over the ob- 
serving run. The other right-hand terms are the phase 
offset 



t'rkn = '/'r"*'''^ + <ArTn of an arbitrarily chosen ref- 
erence antenna r S {«}, and the phase ambiguity term 

€kn^ = "^^^^kn with integer N^^.^ that maps into 
the [0, 27r) domain. 

The antenna-based phase corrections are split into in- 
strumental and ionospheric parts on the basis of their tem- 
poral and spatial behaviour. The phase corrections are 
filtered by iterative estimation of invariant instrumental 
phases (together with the phase ambiguities) and time- and 
space-variant ionospheric phases. The instrumental phases 
are estimated by robust averaging (-1-3 a rejection) over all 
time intervals n: 



Xinstr / i j,cal 

ikn 



mod 2tt\ 



The phase ambiguity estimates follow from 



Tambig 



= 27r round 



/27r 



(12) 



(13) 



where the round() operator rounds a number to the nearest 
integer value. The instrumental phase offset of the reference 
antenna is arbitrarily set to zero. The ionospheric phases 
are constrained by fitting a time-varying spatial gradient 
to the phases over the array. The gradient fit consists 
of an initial estimate directly from the calibration phase 
corrections, followed by a refined fit by using the LM solver 
to minimize 



i:f( 



^kn ' (^i 



(14) 



where aj^ is the position of antenna i. The ionospheric phase 
offset of the reference antenna is arbitrarily set to zero, 
which makes it a pivot point over which the phase gradi- 
ent rotates. Higher order ionospheric effects are assumed to 
average to zero in Equation 1121 

3.2. Initial Phase Calibration and Initial Sky Model 

The instrumental phase calibration method described in 
Section 13.11 assumes that the time-averaged ionospheric 
phase gradient over the array in the direction of the bright 
phase calibrator is zero. Any non-zero average is absorbed 
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into the instrumental phase estimates, causing a posi- 
tion shift of the whole target field and thereby invalidat- 
ing the astrometry. Before entering the calibration cycle 
(Sections l3.3H3.5p . SPAM requires restoration of the as- 
trometry and determination of an initial sky model and 
initial ionospheric calibration. 

To restore the astrometry, the instrumentally corrected 
target field data from Sect ion [3TT] is phase calibrated against 
an apparent sky model (AIPS task CALIB). The default is a 
point source model, using NVSS catalog positions (Condon 
et al. [T994lfT998l) . power-law interpolated fluxes from N VSS 
and WENSS/WISH catalogs (Rengehnk et al. ITM7|) and 
a given primary beam model. To preserve the instrumen- 
tal phase calibration as obtained in Section [XT] during fur- 
ther processing, time-variable phase corrections resulting 
from calibration steps in this and the following sections are 
stored in a table (AIPS SN table) rather than applied di- 
rectly to the visibility data. The sky model calibration is 
followed by wide-field imaging (AIPS task IMAGR) and 
several rounds of phase-only self-calibration (CALIB and 
IMAGR) at the highest possible time resolution, yielding 
the initial sky model and initial phase calibration. 

For wide-field imaging with non-coplanar arrays, the 
standard imaging assumptions that the relevant sky area 
is approximately flat and the third baseline coordinate (w- 
term in Equation [7]) is constant across the FoV are no 
longer valid. To overcome this, SPAM uses the polyhedron 
method (Perley I1989al Cornwell & Perley [T99^ that di- 
vides the large FoV into a hexagonal grid of small, par- 
tially overlapping facets that individually do satisfy the as- 
sumptions above (AIPS task SETFC). Additional facets are 
centered on relatively bright sources inside and outside the 
primary beam area t o reduce image artefacts due to pixel- 
lation fPerlev [T989bl Briggs & CornweU[199 3 Brig gs[T995| 
Voronkov & WieringaHTOl Cotton & Uson [W7|) . 

The Cotton-Schwab algorithm (Schwab 119841 Cotton 
119991 Cornwell, Braun & Briggs I1999P is a variant of 
CLEAN deconvolution (HogbomlMl Clark [HSDl) that al- 
lows for simultaneous deconvolution of multiple facets, us- 
ing a different dirty beam for each facet. Boxes are used to 
restrict CLEANing to real sky emission, making sure that 
sources are deconvolved in the nearest facet only (CLEAN 
model components are stored in facet-based AIPS CC ta- 
bles). After deconvolution, the CLEAN model is restored 
to the relevant residual facets (AIPS task CCRES) using a 
CLEAN beam, and the facets are combined to form a single 
image of the full FoV (AIPS task FLATN). 

3.3. Peeling 

To construct a model of ionospheric phase rotations in ar- 
bitrary viewing directions within the FoV, SPAM requires 
measurements in as many directions as possible. When no 
external sources of ionospheric information are available, 
the target field visibilities themselves need to be utilized. 
Calibration on individual bright sources in the FoV can sup- 
ply the required information, even in the presence of higher 
order phase structure over the array. After instrumental 
phase offsets are removed, phase calibration corrections are 
an relative measure of ionospheric phase: 



cal 
ikn 



'Prkn Vikn 



(15) 



SPAM uses the peeling technique (Noordam I2004p to 
obtain phase corrections in different viewing directions. 
Peeling consists of self-calibration on individual sources, 
yielding per source a set of time-variable antenna-based 
phase corrections and a source model, after which the 
source model is subtracted from the visibility data set while 
temporarily applying the phase corrections (AIPS tasks 
SPLIT, UVSUB and CLINV/SPLIT). 

For peeling to converge, the source needs to be the 
dominant contributor of flux to the visibilities on all base- 
lines. Especially at low frequencies, the presence of many 
other sources in the large FoV adds considerable noise 
to the peeling phase corrections. To suppress this effect, 
the following steps are performed: (i) The best available 
model of the apparent sky is subtracted from the visibil- 
ity data while temporarily applying the associated phase 
calibration(s). The initial best available model and asso- 
ciated phase calibration is the self-calibration output of 
Section 13.21 Individual source models are added back be- 
fore peeling, (ii) Sources are peeled in decreasing flux order 
to suppress the effect of brighter sources on the peeling of 
fainter sources, (iii) Calibration only uses visibilities with 
projected baseline lengths longer than a certain threshold. 
This excludes the high 'noise' in the visibilities near zero- 
length baselines from the coherent flux contribution of im- 
perfectly subtracted sources. 

The radio sky can be approximated by a discrete num- 
ber of isolated, invariant sources of finite angular extend. 
Visibilities in the ME (Equation [7]) for a single integration 
time n can therefore be split into a linear combination of 
contributions from individual sources k: 



(16) 



k k •' •' 

-27rj[Mij„i+t).j-„m+iUij„(n-l)] dm 



The subtraction of all but the peeling source k' from the 
measured visibilities in step (i) above can be described as 



V, 



ijk' n 



k^k' 



t \ 1 T^model 
•. 9jkn) ^ij 



^ i]kn 



(17) 



with g, 



ikn 



9^(1 



ki '"fc' 



Wikn the best available cali- 
bration in the viewing direction of source fc, and V^Jfcrf"' ^^^^ 
visibilities that are derived from the best available model 
^fjk^''^ of source k. The peeling itself consists of iterative 
calibration and imaging steps of the peeling source k' . The 
calibration (Equation [TO] with p = 1) updates the antenna 
gain corrections by minimizing 



EE- 



T/ mo del 
ijn II * ijk' n 



9in 9 jn^ijk' n II i 



(18) 



where we used Equation II II with 



= C^*' = 0. 



while the imaging step updates 1™^^''^ and therefore V^§'lf'. 

In practise, due to incompleteness of the sky model 
and inaccuracies in the phase calibration, there will always 
remain some contaminating source flux in the visibilities 
while peeling. Complemented with system noise, sky noise, 
residual RFI and other possible sources of noise, the noise 
in the visibilities propagates into the phase corrections from 
the peeling process. 
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Absolute astromctry is not conserved during peeling, 
because self-calibration allows antenna-based phase correc- 
tions to vary without constraint. In subsequent peeling cy- 
cles, small non-zero phase gradients in the phase residu- 
als after calibration can cause the source model to wander 
away from its true position. In SPAM, astrometry errors 
are minimized by re-centering the source model to its true 
(catalog) position before calibration in each self-calibration 
loop. By default, SPAM re-centers the peak of the model 
flux to the nearest bright point source position in the NVSS 
catalog (Condon et al. 119941 11998p . It is recommended to 
visually check the final peeling source images for possible 
mismatches with the catalog (e.g., in case of double sources 
or sources with a spatially varying spectral index). 

While peeling; SPAM attempts to calibrate sources on 
the highest possible time resolution, which is the visibil- 
ity time grid. The noise in the resulting phase corrections 
depends on the signal-to-noise ratio (SNR) of the source 
flux in the visibilities. To increase the number of peeling 
sources and limit the phase noise in case of insufficient SNR, 
SPAM is allowed to increase the calibration time interval 
beyond the visibility integration time up to an arbitrary 
limit. Through image plane analysis, SPAM estimates the 
required calibration time-interval per source: 



(19) 



where rif is the required number of integration times in a 
calibration interval, is the total number of integration 
times within the observation, a is the minimum required 
SNR per integration time (a tweakable parameter that sets 
the balance between the SNR and the time resolution of 
the peeling phase corrections), and S and (7^ are the mea- 
sured source peak flux and local background noise level in 
the image. For a fixed upper limit on the calibration time 
interval, an increase in a results in a decrease in the num- 
ber of peeling sources. For < 1. phase corrections are 
determined on the visibility time grid. For > 1, a spline 
is used to resample the phase corrections per antenna in 
time onto the visibility time grid. 

Apart from SNR issues, the number of sources that can 
be peeled is fundamentally limited by the available number 
of independent visibility measurements. When peeling 
sources, self-calibration fits Ng{N^ — 1) phase solutions per 
calibration time interval to the visibility data, where is 
the number of antennas. For self-calibration to converge to 
an unique combination of phase solutions and source model, 
this number needs to be much smaller than the number of 
independent visibility measurements. The maximum of vis- 
ibilities measurements that is available in one calibration 
time interval is given by N^{n^)N^{N^ — l)/2, with the 
number of frequency channels and (rif) the average number 
of visibility integration times in a calibration interval. In the 
ideal case, when we assume that each visibility is an inde- 
pendent measurement, the determination of antenna-based 
phase corrections for all peeling sources is well constrained 
if 



(20) 



The applications presented in this article do satisfy this 
minimal condition (see Section [?]). 

Equation [20] is equivalent to stating that the num- 
ber of degrees-of- freedom (DoF; the difference between the 



number of independent measurements and the number of 
model parameters) should remain a large positive number. 
Correlation between visibilities over frequency and time 
may reduce the number of independent measurements dras- 
tically, thereby also reducing the number of DoFs. The ex- 
act number of DoFs for any data set is hard to quantify. 
When this number becomes too low, the data is 'over-fitted' 
(e.g., Bhatnagar et al. l2008p . which could result in an artifi- 
cial reduction of both the image background noise level and 
source flux that is not represented in the self-calibration 
model ( Wieringa I1992p . Although we have found no evi- 
dence of this effect occuring in the applications presented 
in this article, the SPAM user should be cautious not to 
peel too many sources. In case of a high number of available 
peeling sources, one can choose a subset with a sufficiently 
dense spatial distribution over the FoV (e.g., one source per 
isoplanatic patch; see Section [3.5p . 



3.4. Ionospheric Phase Screen Model 

The phase corrections that are obtained by peeling several 
bright sources in the FoV (Section 13. 3p are only valid for 
ionospheric calibration in a limited patch of sky around 
each source. To correct for ionospheric phase errors over 
the full FoV during wide-field imaging and deconvolution, 
SPAM requires a model that predicts the phase correction 
per antenna for arbitrary viewing directions. 

SPAM constructs a quasi-physical phase screen model 
that attempts to accurately reproduce and interpolate the 
measured ionospheric phase rotations (or more accurately: 
the peeling phase corrections). The phase screen is deter- 
mined independently for each visibility time stamp, there- 
fore we drop the n-subscript in the description below. 
Figure [2] is a schematic overview of the geometry of iono- 
spheric phase modeling in SPAM. The ionosphere is repre- 
sented by a curved phase screen at a fixed height h above 
the Earth's surface, compliant to the WGS84 standard 
(NIMA I1997P . The total phase rotation experienced by a 
ray of radio emission traveling along a LoS through the 
ionosphere is represented by an instantaneous phase rota- 
tion (/)'°"(p, C) on passage through the phase screen that is 
a function of pierce point position p and zenith angle C- For 
a thin layer {Ah ^ S; see Section [2?3| . the dependence of 
0ion Qj-^ ^ represented by a simple airmass function, 

so that 



cos(C) 



(21) 



SPAM uses an angular local longitude/latitude coordi- 
nate system to specify p, relative to the central pierce point 
from array center to field center. For the applications pre- 
sented in this article, the angular distances between pierce 
points over the relevant ionospheric domain are all < 5 de- 
grees, which effectively makes the pierce point vector p a 
2-dimensional cartesian vector. 

The 2-dimensional phase screen 0'°"(p) is defined on a 
set of KL base vectors, generated from the instantaneous 
pierce point configuration {Pj^.} and an assumed power-law 
shape for the phase structure function (Section 12. ip . The 
KL base vector generation and interpolation is based on the 
work by van der Tol & van der Veen (|2007p and is described 
in detail in Appendix [Al The phase screen model requires 
one free parameter per KL base vector. The initial com- 
plete set of KL base vectors is arbitrarily reduced in order 
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Figure 2: Schematic overview of the SPAM thin ionospheric 
phase screen model geometry. For clarity, only two spatial 
dimensions and one calibration time interval are consid- 
ered. In this overview, five ground-based array antennas 
(labelled 1 to 5) observe three calibrator sources (colored 
red/green/blue and labelled A to C) within the FoV. The 
(colored) LoSs from the array towards the sources run par- 
allel for each source and pierce the phase screen at fixed 
height h (colored circles). The LoS from antenna i at Earth 
location towards a peeling source k at local sky position 
s^f. intersects the phase screen at a single pierce point p^^, 
under a zenith angle Q^. For a single LoS from antenna 1 
towards source A, we have indicated how the pierce point 
position Pjj, = P]^^ and zenith angle C^fc = Cia relate to 
the antenna position = Si^ and the local sky position 
s^^ of the source. For some LoSs the pierce points 



may overlap (or nearly overlap), as is the case for IC & 4A 
and 2C & 5A in our example. The total (integrated) phase 
rotation along any LoS through the ionosphere is modeled 
by an instantaneous phase rotation at the phase screen 
height. For example, radio waves traveling along LoSs from 
source A towards antennas 1 to 5 experience an instanta- 
neous phase rotation </)-™ = to (/ig^, respectively, while 
passing the screen at their related pierce points p^j. = pj^^ 
to P5AJ respectively. Peeling the three calibrator sources 
yields measurements of the ionospheric phases relative 
to a common reference antenna (in this example antenna 
3; encircled). 



by selecting a subset based on statistical relevance (princi- 
ple component analysis). This reduces the effect of noise in 
the peeling solutions on the model accuracy and simulta- 
neously limits the number of model parameters. However, 
the subset should still be large enough to accurately re- 
produce the peeling phase corrections. Per visibility time 
stamp, the KL base vectors are stored for later use during 
imaging (for this purpose, we mis- use the AIPS OB table). 
As an example, the first six interpolated KL base vectors 



for a single configuration of ionospheric pierce points are 
plotted in Figure [31 

The peeling phase corrections are interpreted to 
be relative measurements of the absolute ionospheric phase 
screen model (/)'°" (p, Q which may be determined up to a 
constant. The model parameters are determined by mini- 
mizing the differences between the observed and the model 
phases using the LM non-linear least-squares solver, for 
which a sum needs to be defined. From Equation [15] 
it follows that 



ambig 



(22) 



Consequently, the phase correction in the direction of 
source k for a baseline consisting of antennas i and j is 



The sum is defined as: 



-EEE 



i j>i 



/ambig 



(23) 



[r''iP^kX^k) - 0'™(p,fc, C,fc)] mod 2n 



(24) 



This definition has several properties: (i) By remapping the 

terms into the [0, 2tt) domain, the phase ambiguity terms 
do not have to be fitted explicitly, (ii) the terms of all 
calibrator sources are weighted equally, so the model is not 
biased towards the brightest source (as is the case for self- 
calibration), and (iii) using terms from all possible an- 
tenna pairs prevents a bias towards the reference antenna. 

Using Equation [211 the LM solver yields a set of model 
parameters per visibility time stamp. These are stored for 
later use during imaging (AIPS NI table). The square root 
of the average of the x^ terms equals the average RMS 
phase residual between peeling and model phases. Time in- 
tervals that have a bad fit are identified and removed by 
means of an upper limit (-1-2.5 cr rejection) on the distribu- 
tion of RMS phase residuals over time. 

Convergence of the LM solver is troubled by 27r phase 
ambiguities, because these introduce local minima in 
space. A good initial guess of the model parameters greatly 
helps to overcome this problem. To this purpose, SPAM es- 
timates the global phase gradient over all the pierce points 
directly from the phase corrections 0^^' and projects it onto 
the KL base vectors before invoking the LM solver. 

Figure [H shows an example of an ionospheric phase 
screen that was constructed as described above. The pierce 
point layout consists of multiple projections of the ar- 
ray onto the phase screen. The low density of calibra- 
tors causes a minimal overlap between array projections. 
Figure \5\ shows a comparison between time-sequences of 
phase corrections from self-calibration, peeling and model 
fitting. Because the self-calibration corrections are a flux- 
weighted average for the full FoV, they are biased towards 
the brightest source. They look somewhat similar to the 
peeling solutions of the brightest source, but the latter con- 
tains additional fluctuations that vary on a relatively short 
timescale. The model phases appear similar to the peeling 
phases, but vary more smoothly. Their values fall some- 
where in between the self-calibration phases and the peel- 
ing phases. The difference between the peeling phases and 
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Alan^ritiidc [dcg] 

Figure 3: Plots of the interpolations of the first six KL base vectors, derived for an artificial but realistic configuration 
of ionospheric pierce points. In this example, the pierce points (black crosses) are calculated for a single time instance 
during a 74 MHz VLA-B observation with 13 available calibrator sources in the ~ 10 degree FoV, adopting a phase 
screen height h = 200 km and a structure function power-law slope 7 = 5/3. The horizontal and vertical axes represent 
angular distances in East- West and North-South directions, respectively, as seen from the center of the Earth, relative 
to the phase screen's pierce point along the line-of-sight from array center to pointing center, with East- and Northward 
offsets being positive. At this height, a 0.1 degree angular offset represents a physical horizontal offset of ^ 11.5 km. 
The direction-dependent phase for each interpolated KL base vector is color-coded and scaled to an arbitrary amplitude 
range. 



model phases are mainly caused by the constraints on the 
spatial variability of the phase screen model. 

3.5. Imaging 

With an ionospheric phase screen model available for 
a given visibility data set, antenna-based phase correc- 
tions for any direction in the wide FoV can be calculated 
(Equation [22) . Because each visibility consist of contribu- 
tions from visible sources in different viewing directions, 
there is no simple operation that removes the ionospheric 
phase rotations from a visibility data set prior to imag- 
ing. Instead, SPAM requires an algorithm that calculates 
and applies the appropriate model phase corrections during 
imaging and deconvolving for different parts of the FoV. 

SPAM works under the assumption that there exists a 
fixed angular isoplanatic patch size on the sky, with a pro- 
jected size at ionospheric height smaller than the scale size 
of ionospheric phase fluctuations, over which variations in 
ionospheric phase rotation are negligible. Each isoplanatic 
patch requires at least one phase correction per antenna per 
visibility time interval. For the VLA at 74 MHz, the isopla- 



natic patch size is estimated to be 2-4 degrees (Cotton & 

Condon unnn). 

The facet-based polyhedron method for wide-field imag- 
ing (see Section 13. 2p allows for a relatively simple imple- 
mentation of ionospheric phase correction (Schwab I1984p . 
By choosing a facet size smaller than the isoplanatic patch 
size, a set of model phase corrections calculated for the 
center of a facet are assumed to be accurate for the whole 
facet area. Ionospheric phase model corrections are calcu- 
lated and stored (AIPS SN tables) for each facet center in 
the FoV prior to imaging and deconvolution. For the ad- 
ditional facets centered on bright sources (see Section [3T2)) . 
model phase corrections are optionally replaced by peeling 
phase corrections to allow for optimized calibration towards 
these sources. 

The SPAM imaging and deconvolution procedure is 
similar to the procedure used for the field-based calibra- 
tion method by Cotton et al. (|2004p . which differs from 
the standard Cotton-Schwab algorithm by the temporary 
application of the facet-based phase corrections (AIPS 
tasks SPLIT and CLINV/SPLIT) to the visibility data for 
the duration of major CLEAN cycles on individual facets 
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Figure 4: Example of an ionospheric phase screen model fit. 
The color map represents an ionospheric phase screen at 
200 km height that was fitted to the peeling phase solutions 
of 8 calibrator sources at time-interval = 206 of 10 sec- 
onds during a VLSS observing nm of the 74 MHz VLA in 
BnA-configuration (see Section[H the J1300-208 data set). 
The plot layout is similar to Figure [3l The overall phase 
gradient (depicted in the bottom-left corner) was removed 
to make the higher order terms more clearly visible. The 
collection of pierce points from all array antennas to all 
peeling sources are depicted as small circles., The color in 
the circle represents the measured peeling phase (the refer- 
ence antenna VLA N36 was set to match the phase screen 
value). The size of the circle scales with the magnitude of 
the estimated phase residual after model correction. The 
overall RMS phase residual (7^^^^^^ = 21.799 degrees (av- 
eraged over all pierce points) was one of the better fitting 
results during this particular observing run. 



(AIPS tasks IMAGR and UVSUB). After deconvolution, 
facets are combined to form a single image of the full FoV 
(AIPS task FLATN). Because antenna-based phase correc- 
tions change very little between adjacent facets, the com- 
plete set of partly overlapping facet images combine into a 
continuous image of the FoV. 



4. Applications 

To demonstrate the capabilities of SPAM, we have defined 
three test cases based on observations with the VLA at 
74 MHz (Kassim et al. \^UU7\\ . In each test case, SPAM 
is used for ionospheric phase calibration and imaging of a 
VLSS visibihty data set (Cohen et al. l2007p . following the 
steps described in Section [31 In the first test case, SPAM 
was applied to simulated data to validate basic functional- 
ity in a controlled environment. In the next two test cases, 
SPAM was applied to visibility data from real observations 
under varying ionospheric conditions. We compare SPAM 
performance against self-calibration (SC) and field-based 
calibration (FBC) by analyzing the resulting images. The 
setup and results of these test cases are described in detail 
in the following sections. 
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Figure 5: Example of phase corrections from different steps 
in the ionospheric calibration process, resulting from pro- 
cessing a VLSS data set with SPAM (see Section [H the 
real J0900-I-398 data set). The antenna under considera- 
tion is VLA E28, with W20 being the reference antenna 
(an 5.7 km east-west baseline). The plots represent 25 min- 
utes of observing time, using a 10 second time resolution. 
Top: Antenna-based phase corrections resulting from self- 
calibration on the whole FoV. Middle: Phase corrections re- 
sulting from peeling the brightest (30 Jy) source. Bottom: 
Corrections resulting from ionospheric phase modeling in 
the direction of the (same) brightest source. 



4.1. Data Selection, Preparation and Processing 

In this Section we describe how the visibility data 
sets for the three test cases were selected/constructed. 
Furthermore, we present details on how these data sets were 
processed by SPAM into calibrated images of the FoV. 

Two VLSS observations, at pointing centers J0900-I-398 
and J1300-208, respectively, have been picked from more 
than 500 available VLSS observations on the following cri- 
teria: (i) both fields contain a relatively large number of 
bright sources that can serve as calibrators, and (ii) the 
ionospheric conditions during the observations appear to 
be relatively good (J0900-h398) and relatively bad (J1300- 
208). The presence of more than 5 bright sources of at 
least 5 Jy compensates for the relatively poor efficiency 
of the VLA 74 MHz receiving system (Kassim et al. I2007|) . 
The ionospheric conditions were derived from the appar- 
ent smearing of point sources in the images, due to resid- 
ual phase errors after applying FBC. From experience, we 
adopted the qualification 'good' when the mean width of 
apparent point sources was at most 5" larger than the in- 
trinsic 80" resolution, while for 'bad' conditions the mean 
point source width was larger by at least 15". In terms of 
Strehl ratio R (Equation [B]); 'good' and 'bad' conditions 
correspond with R > 0.996 and R < 0.966, respectively. 
Additionally, candidate fields were visually inspected for 



H.T. Interna et al.: Ionospheric Calibration of LF Radio Observations 1. 



13 



Field name 


VLSS J0900+398 (simulated) 


VLSS J0900+398 (real) 


VLSS J1300-208 (real) 


Pixel size*^ 


18.9" 


18.9" 


11.1" 


Number of facets 


347 


243 


576 


Facet separation 


1.18° 


1.18° 


0.62° 


SPAM calibration cycles'' 


1 


1 


3 


Peeling sources 


10^ 


11 


9 


KL model height 


1000 km"* 


200 km 


200 km 


Fitted KL model terms 


15 


15 


20*= 


Rejected time intervals 


/ 464 


25 / 464 


86 / 484 


Model fit phase RMS 


3.0° ±0.8° 


21.3° ±2.4° 


23.2° ±3.2° 


Peeling corrections applied directly 


no 


yes 


yes 



" The pixel size for all field-based calibration images is 20" . 

'' Adding more cycles did not significantly improve the image quality. 

Arbitrarily limited to mimic a more realistic scenario. 

Increased to improve match with FBC phase screen. 
" In this case, 15 terms proved to be insufficient. 



Table 1: Overview of processing parameters for the three data sets that are handled with SPAM as defined in the test cases. 



evidence of residual phase errors by the presence or ab- 
sence of image artefacts near bright sources, which lead to 
the final selection of the two fields mentioned above. 

The difference in observed ionospheric conditions be- 
tween the two real data sets may be the result of the dif- 
ference in array size and elevation of the target field. From 
the VLA site at ±34 degrees declination, the J0900±398 
field was observed in B-configuration (up to 11 km base- 
lines) at relatively high elevation, while the J1300-208 field 
was observed in BnA-configuration (up to 23 km baselines) 
at relatively low elevation. For the J1300-208 observation, 
the array observed through the ionosphere at larger separa- 
tions and along longer path lengths than for the J0900±398 
observation, which is expected to result in both larger and 
less coherent phase errors over the array. 

Because both real data sets have been previously cal- 
ibrated and imaged with FBC, the data sets were al- 
ready partly reduced at the start of SPAM processing. 
Instrumental calibration was applied (including instrumen- 
tal phase calibration, similar to Section 13. ip . most RFI- 
contaminated data was flagged and the spectral resolution 
was reduced (see Cohen et al. l2007l for details), but no FBC 
has been applied yet. For the simulated data set, which is 
based on the real J0900±398 observations, the measured 
visibilities were replaced by noiseless model visibilities of 
an idealized sky, consisting of 91 bright point sources with 
peak fluxes (larger than 1 Jy) and positions as measured 
in the J0900±398 FBC image. For each point source, the 
corresponding model visibility phases were corrupted using 
the direction-dependent ionospheric phase model that was 
obtained with FBC to correct the real J0900±398 data. 

FBC images of the two real data sets were available in 
the VLSS archive. For the simulated J0900±398 data set, 
an 'undisturbed' image was made before applying the iono- 
spheric phase corruptions. All three VLSS data sets have 
been processed with SPAM, yielding both an SC image and 
an ionosphere-corrected SPAM image. Relevant details on 
the processing can be found in Table [TJ For SC and SPAM 
imaging, we adopted most of the imaging-specific settings 
from FBC (like uniform weighting). Noticeable differences 
are the use of CLEAN boxes, a smaller pixel size and a 
different facet configuration. 



By choosing a minimum SNR per time interval of 15 
and a maximum peeling time interval of 4 minutes (see 
Equation [T9|) . SPAM was able to peel 10 sources in each 
of the real data sets. Lowering the SNR resulted in a much 
larger scatter in the peeling phases over time, or prevented 
peeling to converge at all. The peeling time upper limit was 
chosen to roughly match the spatial density of calibrator 
sources used in FBC. Determining phase corrections on a 
4 minute time scale could result in undersampling the time 
evolution of ionospheric phase errors. Note that this only 
applies to the faintest of the calibrator sources. The limi- 
tations on spatial and temporal sampling of the ionosphere 
are dictated by the given sensitivity of the VLA. 

Because of the high SNR, all 91 sources in the simulated 
J0900±398 data set qualified for peeling at the highest time 
resolution of 10 seconds. To mimic a more realistic scenario 
for further SPAM processing, the number of calibrators was 
arbitrarily limited to 10. Generally, for all data sets, the 
images of peeling sources showed larger peak fluxes and 
less background structure than their counterparts in the 
SC image, although the contrast became less apparent for 
weaker and extended (mostly doubles) peeling sources. 

As stated in Section !??^ the number of peeling sources is 
fundamentally limited by the requirement for a large posi- 
tive number of degrees-of-freedom in the available visibility 
data. The minimal requirement is given in Equation 1201 
Typically, for the VLSS data sets, there were 25 active an- 
tennas, 12 frequency channels and 6 visibility intervals (of 
10 seconds) in an average peeling interval of 1 minute. In 
our test cases, we typically peel 10 sources, which is much 
less than 25 x 12 x 6/2 = 900, thereby satisfying the minimal 
requirement. 

Due to the uncertainty in their optimal values, it is left 
to the SPAM user to specify the phase screen model order 
(the number of KL base vectors), the height h of the phase 
screen and the power-law exponent 7 of the phase structure 
function. For the applications presented here, we used h = 
200 km and 7 = 5/3, which is compliant to the measured 
values given in Section 12.11 given the uncertainty in these 
values. For the simulated data set, we chose instead h = 
1000 km to better match the corrupting FBC ionospheric 
phase model that is attached to the sky plane at infinite 
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Figure 6: The grayscale map represents the residual phase 
RMS between the distorting and correcting ionospheric 
phase models across the primary beam area, averaged over 
baselines and time. The phase RMS was calculated for 
a hexagonal grid of viewing directions across the FoV. 
Each viewing direction is depicted by a small circular area. 
Overplotted is a contour map of the point sources as seen in 
the SPAM image (which extends slightly beyond the grid of 
circles). The 10 peeling sources are marked by circles. The 
correspondence between the models is largest near the cal- 
ibrator sources and over a large part of the inner primary 
beam. The discrepancy is largest near the South-East and 
North- West borders, away from the calibrators. 



height. These values gave satisfactory results for the test 
applications presented here, but can be further optimized. 
The optimal model order was found to lie in the range of 
15-20 terms, which is 1.5-2 times the number of available 
peeling sources. Increasing or decreasing the model order 
caused the model fit to be less accurate or more problematic 
in terms of convergence. 

For both the simulated and real J0900-I-398 data sets, no 
improvement in background noise was observed by adding 
a second calibration cycle after the first. This indicates fast 
convergence of the SPAM calibration method for quiet iono- 
spheric conditions, where the initial self-calibration is al- 
ready close to the best achievable calibration of SPAM. For 
the real J 1300-208 data set, adding up to third calibration 
cycle did improve over the previous cycles. 

4.2. Phase Calibration Accuracy 

For the simulated J0900-f 398 data set, the absolute accu- 
racy of ionospheric calibration can be determined by a di- 
rect comparison between the corrupting FBC phase screen 



and the correcting SPAM phase screen. To this purpose, 
phase corruptions and corrections were calculated from the 
models for a hexagonal grid of 342 viewing directions within 
the FoV. Per viewing direction, the RMS phase error was 
calculated by differencing of the phases from both mod- 
els and averaging over all time stamps and baselines. The 
result is depicted in Figure [6l 

For areas near the calibrators and in the center of the 
field in general, there is a relatively good match between 
the input and output model, with typical RMS phase er- 
rors < 5 degrees. The absence of calibrator sources south- 
west of the field center still results in relatively accurate 
predictions by the SPAM model. In the direction of peeling 
sources, the measured RMS phase error can be split into a 
contribution from inaccuracies in the peeling process and a 
contribution from imperfect model fitting. The latter is ap- 
proximately 3 degrees (Table [1]) , therefore the RMS phase 
error introduced by peeling is < 4 degrees. Considering the 
model setup, the only possible source of error is contami- 
nation from other sources while peeling (which appears to 
happen despite the initial subtraction of the SC model). 

Overall, the change in model base from the corrupt- 
ing FBC model (5 Zernike polynomials) to the correcting 
SPAM model (15 KL vectors) has a constant accuracy over 
large parts of the FoV. Towards some parts of the edge of 
the field the phase errors are substantially larger, up to 20- 
25 degrees at worst. This agrees with the different asymp- 
totic behaviour towards large radii of the Zernike model 
(diverge to infinity) and the KL model (converge to zero) 
in the absence of calibrators. The presence of calibrator 
sources near the edge (like the one on the North-East edge 
of the field) leads to a better local match between corrupt- 
ing ionosphere and correcting model. 

For the real observations, in the absence of external 
sources of information (e.g., GPS measurements), it is not 
possible to derive the absolute accuracy of ionospheric cal- 
ibration from the observations themselves. Instead, the 
residual RMS phase error of the model fit to the peeling 
phases is used as an relative indicator for calibration accu- 
racy over time. For both the real J0900+398 and J1300-208 
data sets, the residual RMS phase error of ^ 22 degrees is 
much larger than for the simulated data. This already ex- 
cludes rejected time stamps with exceptionally large RMS 
values. By inspecting model fits on individual time stamps, 
we found that there are often a few pierce point phases 
that deviate significantly more from the fitted model than 
most neighbouring points. These errors do not appear to 
be antenna-based instrumental errors, because peeling solu- 
tions for the same antenna towards other calibrator sources 
do not deviate in the same manner. Typically, these deviat- 
ing points persist for a few time stamps before disappear- 
ing. The ionosphere may be responsible for these very small 
scale deviations. Another possibility is that the peeling so- 
lutions are (sometimes) noisy due to limitations in source 
SNR. 

4.3. Background Noise 

In this and the next sections, we revert to analyzing image 
properties for an indirect, relative comparison between the 
different calibration techniques. In the presence of residual 
phase errors, part of the image background noise level con- 
sists of residual sidelobes after CLEANing. The local side- 
lobe noise increases with both the RMS phase error and the 
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Field name VLSS J0900+398 (simulated) VLSS J0900+398 (real) VLSS J1300-208 (real) 



Mean background noise 


a [mjy beam ^ ] : 






Undisturbed 


3.0 


- 


- 


SC 


10.2 


71 


92 


FBC 


- 


87 


118 


SPAM 


6.7 


67 


68 


Number of sources with 


a peak flux larger than 5ct: 






Undisturbed 


91 


- 


- 


SC 


91 


393 


374 


FBC 


— 


310 


285 


SPAM 


91 


372 


392 


5(7 source fraction with 


an NVSS counterpart within 80" 






Undisturbed 


1. 






SC 


1. 


0.83 


0.60 


FBC 




0.86 


0.74 


SPAM 


1. 


0.97 


0.97 



Table 2: Overview of results from calibrating and imaging three test case data sets with no ionosphere 
(Undisturbed), self-cahbration (SC), field-based cahbration (FBC) and SPAM. 



local source flux. When measured over a large image area, 
the mean sidelobe noise depends mainly on RMS phase er- 
ror. For all relevant output images, the mean image noise 
cr was determined by fitting a Gaussian to the histogram 
of image pixel values from the inner quarter radius of the 
FoV (AIPS task IMEAN). Note that these images have not 
been corrected for primary beam attenuation. The results 
are given in Table [2] 

Because no noise was added to the simulated J0900+398 
data set, the resulting image noise of 3.0 mJy beam^^ in the 
undisturbed image is caused by incomplete UV coverage 
and inaccuracies in the imaging process (see Section l3.2p . 
limiting the dynamic range to ^ IQ'^. The local noise is 
highest near the sources, but significantly less near the 
brightest 10 sources with dedicated facets centered on their 
peak position. The SC and SPAM images from this data 
set were created using the same facet configuration. The 
SC image noise of 10.2 mJybeam~^is 3.4 times as high 
as the undisturbed image noise, therefore dominated by 
phase error induced sidelobe noise. The SPAM image noise 
of 6.7 mJybeam~^is a significant improvement over the 
SC image, but still 2.2 times as high as in the undisturbed 
image. The local noise in the SC and SPAM images has 
increased most apparently near bright sources as compared 
to the undisturbed image, which confirms the presence of 
residual phase errors after calibration. 

For the real J0900-h398 data set, both the SC and SPAM 
images have an image noise of ^ 70 mJybeam"^ . The 
SPAM image noise is slightly lower than SC. The local noise 
in the SC image is higher near bright sources. This is not 
the case in the SPAM image, which must be a direct result 
of an improved calibration accuracy near these sources. The 
FBC image noise for this data set is ~ 20 percent higher, 
a combination of a higher average noise over the FoV and 
higher local noise near bright sources. 

For the real J1300-208 data set, the SPAM image has 
the same image noise as for the real J0900-t-398 data set, 
with no apparent increase near bright sources. At the same 
time, the noise levels in the SC and FBC images have in- 
creased with 30 and 35 percent, respectively. The noise in 



the SC image is highest near the bright sources. The FBC 
noise is highest near the brightest source and remains high 
in the rest of the image. The significant increase of the 
average FBC noise level indicates a dependence on iono- 
spheric conditions, and therefore on calibration accuracy. 
The SPAM image noise appears to have little or no depen- 
dence on varying ionospheric conditions (Figure [7|). 

4.4. Source Properties 

The presence of residual phase errors changes the apparent 
distribution of flux of a source (see Section 12. 2p . In the 
time-averaged image, sources may appear offset from their 
intrinsic position, may suffer from smearing or deformation, 
and sidelobes may be misidentified as sources. Comparing 
the properties of the same sources in differently calibrated 
images allows for a relative comparison of the performance 
of the different calibration techniques. 

To allow for comparison of source properties, we ap- 
plied the source extraction tool BDSM (Mohan et al. I2008P 
on all relevant images. BDSM performed a multiple 2- 
dimensional Gaussian fit on islands of adjacent pixels with 
amplitudes above a specified threshold based on the lo- 
cal image noise cr^ in the image. Multiple overlapping 
Gaussians were grouped together into single sources. We 
applied BDSM to all images, using the default extraction 
criteria, except for the following: a source detection requires 
at least 4 adjacent pixel values above 2.5 (Tlj with at least 
one pixel value above 4(Tl. 

4.4.1. Source Counts 

Due to the non-Gaussian character of the phase-induced 
sidelobe noise, the source catalogs will contain spurious de- 
tections. To suppress these, we removed sources with a peak 
flux smaller than 5 a from the catalogs. The remaining num- 
ber of catalog entries are listed in Table [2] Additionally, 
each catalog was cross-associated against the NVSS cat- 
alog, which has a slightly higher resolution (45"). For an 
average spectral index of —0.8, the NVSS detection limit is 
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Figure?: Grcyscalc plots of a 3.5 x 3.5 square degree area in the VLSS J1300-208 field centered on the bright (40 Jy) 
point source 3C 283. All three images have contours (black lines) overplotted at [0.15, 0.48, 0.83, 1.16, 1.50] Jy. Left: 
Image after self-calibration, middle: image after field-based calibration, and right: image after SPAM calibration. 
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Figure 8: Peak flux ratios of point sources in the simulated 
J0900-I-398 field. Peak fluxes were measured in the self- 
calibration image and the SPAM image, corrected for a 
small CLEAN bias and divided by the input model peak 
fluxes. The size of each dot scales with the input model peak 
flux, ranging from 1.02 to 26.7 Jy. Ideally (without phase 
errors), the peak flux ratios would be scattered around one 
(solid lines) due to image noise dependent errors in the peak 
flux determination. Instead, the peak flux ratio distribu- 
tions along the x- and y-axis are centered around 0.995 and 
0.996, respectively (dotted lines), which is a direct result of 
the residual phase errors. The smaller and larger scatter in 
distribution of SPAM- and self-calibration peak flux ratios 
is consistent with peak flux determination inaccuracies due 
image background noise levels. 



at least 10 times lower than for the VLSS. At the risk of 
missing an incidental ultra-steep spectrum source, we de- 
termined the source fraction that has an NVSS counterpart 
within an 80" radius (one VLSS beamsize), which are also 
listed in Table [21 

For the simulated J0900-I-398 data set, all 91 input 
sources are detected and matched against NVSS counter- 
parts, regardless of the calibration method. Due to the low 
noise levels and the lower limit of 1 Jy on the input source 
catalog, all sources are effectively > 100 a detections. None 
of the sources had more than one Gaussian fitted to it, 
despite the freedom to do so. 

For the real J0900-f398 data set, the higher cr in the 
FBC image is reflected in a smaller number of source de- 
tections as compared to SC and SPAM. SC detects shghtly 
more sources than SPAM, despite the slightly higher a. 
However, there is a very large fraction of sources in the 
SPAM catalog that has an NVSS counterpart, significantly 
larger than for both the SC and FBC catalogs. This sug- 
gests that the SPAM catalog is much less contaminated by 
false detections than the SC and FBC catalogs, resulting 
in a larger absolute number of true detections. 

This is further strenghtened by the results from the real 
J1300-208 data set. For this test case, the SPAM image has 
the largest number of source detections. Again, the SPAM 
catalog has the largest fraction of associations with the 
NVSS catalog, the same fraction as with the J0900+398 
data set. In contrast, the fraction of NVSS counterparts 
for SC and VLSS have both gone down. This is best ex- 
plained by an increase in (non-Gaussian) sidelobe noise in 
the image background due to calibration errors, which cor- 
responds with the observed increase in a . 

4.4.2. Source Peak Fluxes 

The presence of residual phase errors after calibration 
can cause an unresolved source shape to deviate from a 
point source shape. The source flux is redistributed over a 
larger area and the peak flux of the source drops. At 80" 
resolution, most sources in a VLSS field are unresolved. 
Therefore, a mean increase of source widths over the point 
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Figure 9: Peak flux ratios in the simulated J0900+398 field: Left: Peak flux ratios of the 91 extracted sources from the 
SPAM image as compared to the input model sources, plotted as a function of the residual RMS phase error after SPAM 
calibration. Overplotted is the theoretical Strehl ratio (solid line) as given in Equation [HI For larger RMS phase errors, 
the measured peak flux ratios do not follow the theoretical strehl ratio curve. This indicates that systematic phase errors 
dominate the larger RMS phase errors. Right: Same peak flux ratios plotted as a function of absolute position offset 
between extracted sources in the SPAM image and the input model (see Figure [T^ . The presence of a strong correlation 
indicates that residual phase gradients dominate the larger RMS phase errors. 



source width is a direct measure of ionospheric conditions. 
This argument was used in the pre-selection of data sets for 
our test cases. 

For significant source deformations or low SNR sources, 
determination of the shape of individual sources is subject 
to large uncertainties (e.g., Condon et al. I1997p . Because 
determination of peak fluxes is much more robust, we use 
these for a relative comparison of calibration accuracy. 
Starting with the original catalogs as produced by BDSM, 
we associate sources between the undisturbed, FBC, SC 
and SPAM catalogs that lie within 80" of the same NVSS 
source and has a peak flux larger than 5 cr in at least one 
of the two catalogs. 

For the simulated J0900+398 data set, the true peak 
fluxes of all 91 sources are known. A comparison between 
peak fluxes from the undisturbed image and the input 
catalog identifies a small (< 1 percent) CLEAN bias of 
3.6 mJybeam-i ( e.g C ondon et al. [TOMl HMEl Becker, 
White & Helfand I1995P . Ignoring the image noise depen- 
dency of CLEAN bias, we applied this small correction to 
the peak fluxes in the undisturbed, SC and SPAM source 
catalogs before proceeding. Figure [5] shows a comparison of 
the measured-to-input peak flux ratios for sources in the SC 
and SPAM images. The mean peak flux ratio for both im- 
ages is approximately equal and just slightly smaller than 
one. The larger scatter in the SC peak fluxes is consistent 
with a higher a. Using Equation [HI the random part of 
the mean RMS phase error for both SC and SPAM is es- 
timated at 5-6 degrees. This value is comparable to the 
observed RMS phase error over large parts of the SPAM 
image (Section 14. 2p . 

To study the nature of residual RMS phase errors after 
application of SPAM, we plot the RMS phase errors at the 
source positions from Figure [6] against SPAM-to-input peak 
flux ratios (Figure [5]). For Gaussian random phase errors, 



the peak flux ratio is expected to decrease with increased 
RMS phase error as described in Equation [6l However, the 
discrepancy between the data points and Equation [5] in- 
dicates that for larger RMS values the phase errors are 
predominantly systematic rather than random. 

For the real J09004-398 data set, FigurefTUlshows a com- 
parison of peak fluxes for associated sources in the SC, FBC 
and SPAM catalogs. There is a good match between peak 
fluxes measured in the SC and SPAM catalogs. For high 
SNR sources with a peak flux above 1 Jy, the SPAM peak 
fluxes match on average within 1 percent with the SC peak 
fluxes. Similarly, SC and SPAM peak fluxes are on average 
10 percent higher than FBC peak fluxes. The systematic 
increase of peak fluxes for SC and SPAM as compared to 
FBC for many more than the calibrator sources denotes 
a more accurate calibration over large parts of the FoV. 
Towards the low flux end, source detections are slightly bi- 
ased towards the image with the highest noise level, which 
is the FBC image. 

Figure [TT] shows the same comparison of peak fluxes for 
the real J1300-208 data set. For high SNR sources with a 
peak flux above 1 Jy, the SC peak fluxes are by far the 
smallest, while FBC and SPAM peak fluxes are on aver- 
age higher by 15 and 24 percent, respectively. The relative 
loss of peak flux in the SC image is a clear indication of 
the break-down of the assumption of isoplanaticity across 
the FoV. Under the conditions that clearly need direction- 
dependent corrections, the SPAM peak fluxes are on aver- 
age 7 percent higher than the FBC peak fluxes. 

4.4.3. Astrometry 

When the time-average of residual phase errors towards a 
source contains a non-zero spatial gradient, the source will 
appear to have shifted its position in the final image (see 
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FigurclO: Peak fluxes in tfie real J0900+398 field Left: Peak fiux comparison for 367 sources detected in both the self- 
calibration and SPAM images. The straight diagonal line represents equality, the dashed lines represent Bctq deviations 
(where aQ is the combined noise level from both images), and the dotted lines indicate the 5 a detection limit. For bright 
sources (peak fluxes > 1 Jybeam"^ ), the average peak flux ratio is 1.00. Middle: Same for 329 sources in the field-based 
calibration (VLSS) and SPAM images. The average bright peak flux ratio of SPAM over field-based calibration is 1.10. 
Right: Same for 313 sources in the self-calibration and field-based calibration (VLSS) images. The average bright peak 
flux ratio of self-calibration over field-based calibration is 1.10. In all plots, the image noise causes a larger scatter in the 
peak flux determinations of faint sources (< 1 Jybeam^^ ) and consequently, a selection bias towards positively enhanced 
peak fluxes that increases with image noise. 




Figure 11: Peak fluxes in the (real) J1300-208 fleld: Left: Peak flux comparison for 247 sources detected in both the 
self-calibration and SPAM images. For bright sources (peak fluxes > 1 Jybeam"^ ), the average peak flux ratio of SPAM 
over SC is 1.24. Middle: Same for 278 sources in the fleld-based calibration (VLSS) and SPAM images. The average 
bright peak flux ratio of SPAM over fleld-based calibration is 1.07. Right: Same for 202 sources in the self-calibration and 
fleld-based calibration (VLSS) images. The average bright peak flux ratio of field-based calibration over self-calibration 
is 1.15. 



Section [2.2p . This gradient may indicate a limitation of the 
calibration model to reproduce the ionospheric phase cor- 
ruptions (e.g., in the absence of nearby calibrators), but 
may also be introduced by the peeling process. The latter 
occurs when a peeling source is re-centered to the wrong 
catalog position (see Section I3.3P . Because such an error 
propagates into the calibration model, many sources in the 
vicinity of the peeling source may also suffer from a sys- 
tematic astrometric error. 

For the simulated data set, the peak positions of sources 
as determined by BDSM were compared against the posi- 
tions of counterparts in the input model. For the real data 



sets, we compared against the NVSS catalog instead. When 
comparing against NVSS positions, apparently large po- 
sition offsets may occur due to resolution differences and 
spectral variation across the source. Averaged over a large 
number of sources, these offsets should have no preferen- 
tial orientation. In contrast, a residual phase gradient in a 
certain viewing direction is expected to cause systematic 
offsets for groups of sources in a certain preferential direc- 
tion. 

For the simulated J0900-I-398 data set. Figure fT2l shows 
that the positions for both SC and SPAM are accurate to 
within ~ 10", except for a small tail of ^ 15 SPAM sources 
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Figure 12: Position offsets in the simulated J0900+398 field: Left: Offsets between the measured source positions in the 
self-calibration image as compared to the input model. Right: Same for the SPAM image. In both cases, the distribution 
around the origin is non-Gaussian. For the SPAM image, the tail of points extending roughly northwards indicates the 
presence of persistent phase gradients in local parts of the SPAM image. All source position offsets fall well within the 
size of the 80" restoring beam (dotted line). 



that have somewhat larger offsets. These sources are all 
positioned near the edge of the FoV, where the RMS phase 
error is large (Figure [6|) . Figure [9] also confirms this by the 
clear correlation between RMS phase error and absolute 
position offsets. 

For the real J0900+398 data set, the source position 
offsets for SC, FBC and SPAM relative to NVSS cata- 
log positions are plotted in Figure [TSj The larger scatter 
as compared to the simulated J0900+398 data set can be 
the (combined) result of less accurate position measure- 
ments due to higher image noise, resolution and spectral 
differences between the observations and the NVSS catalog 
or larger residual RMS phase errors after calibration. The 
observed scatter for SC is centered around a point that 
is offset from the origin by ~ 5", which is either caused 
by inaccuracies in the initial sky model or during the self- 
calibration process fSection lS.'ip . The scatter of both FBC 
and SPAM offsets is centered close to the origin. The RMS 
of the scatter around the mean position offset is 10.5" for 
both FBC and SPAM (despite the apparently larger scatter 
for SPAM, which is due to a larger number of data points), 
both smaller than the 11.9" for SC. 

For the real J1300-208 data set. the source position off- 
sets for SC, FBC and SPAM relative to NVSS catalog po- 
sitions are plotted in Figure [141 The position scatter for 
all three methods is significantly larger than for the real 
J0900+398 data set, and all suffer from systematic position 
offsets in varying degrees of severity. The position offsets in 
the SC image have a seriously distorted distribution, which 
includes a large tail of points that extends roughly south- 
wards. This indicates the presence of varying systematic 
source offsets over the whole FoV. The distribution of po- 
sition offsets in the FBC image is more compact but also 
asymmetric, and is approximately centered around a point 
that is ~ 10" offset in northward direction from the origin. 
A large number of the SPAM position offsets are clustered 
near the origin, similar to the real J0900+398 data set, but 
there is an additional tail of points that runs roughly north- 



wards. The RMS of the scatter around the mean position 
offset is 20.7", 16.5" and 14.8" for SC, FBC and SPAM, 
respectively, which confirms the apparently strongest clus- 
tering of points in the SPAM position offset plot. 

Systematic position offsets in the images can be reduced 
by distortion and regridding of the images. To this purpose, 
Cohen et al. (|2007p fit a fourth order Zernike polynomial to 
the (time constant) position offsets of typically more than 
100 sources in the FBC images of the VLSS. They estimate 
that, after correction, the final residual position error in the 
full VLSS catalog due to the ionosphere is < 3" in both RA 
and DEC. 

5. Discussion and Conclusions 

The SPAM method for ionospheric calibration has been 
succcsfully tested on one simulated and two carefully se- 
lected visibility data sets of 74 MHz observations with the 
VLA (taken from the VLSS; Cohen et al. [Wf|) . From the 
results of these test cases, we draw the following conclu- 
sions: 

(i) A proof-of-concept is given for several different tech- 
niques that were incorporated in SPAM calibration. The 
peeling technique (Noordam I2004p was succesful in pro- 
viding relative measurements of ionospheric phase errors 
in the direction of several bright sources in the FoV. The 
Karhunen-Loeve phase screen (van der Tol & van der Veen 
I2007|) at fixed height was able to combine these measure- 
ments into a consistent model per time stamp. For rela- 
tively bad ionospheric conditions, it was demonstrated that 
the ionospheric calibration cycle (repeated ionospheric cali- 
bration and subsequent imaging; Noordam l2004p converges 
within a few iterations to a calibration of similar accuracy 
as under relatively good ionospheric conditions (for which 
one iteration was sufficient). 

(ii) Ionospheric calibration with SPAM is more accu- 
rate than the existing self-calibration (e.g., Pearson & 
Readhead I1984p and field-based calibration (Cotton et 



20 H.T. Interna et al.: Ionospheric Calibration of LF Radio Observations I. 



si^lf-calibration position oftsots VLSS position offsets SPAM position offsots 




80 60 40 20 -20 —10 -60 -80 80 60 40 20 -20 -10 -60 -80 80 60 40 20 -20 -10 -60 -80 



ARA |arcsei-| ARA |aiTs«n| ARA |arrsei-| 

Figure 13: Position offsets in the real J0900+398 field Left: Offsets between the measured source positions in the self- 
calibration image as compared to the NVSS catalog. Middle: Same for the field-based calibration (VLSS) image. Right: 
Same for the SPAM image. 
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Figure 14: Position offsets in the real J1300-208 field: Left: Offsets between the measured source positions in the self- 
calibration image as compared to the NVSS catalog. Middle: Same for the field-based cahbration (VLSS) image. Right: 
Same for the SPAM image. 



al. l2004p techniques. Even for relatively compact array con- 
figuration like VLA-B and BnA, significant improvements 
in image quality are obtained by allowing for higher-order 
(i.e., more than a gradient) spatial phase corrections over 
the array in any viewing direction. In the resulting images, 
we obtained dynamic range improvements of 5-45% and 
70-80% under relatively good and bad ionospheric condi- 
tions, respectively. 

(iii) Although the mean astrometric accuracy of source po- 
sitions in SPAM images is similar to or better than for 
self-calibration and field-based calibration, systematically 
larger astrometric errors are present in regions of the out- 
put images of all calibration methods. This is caused by a 
shortage of available calibrators in these regions and posi- 
tional inaccuracies in the reference source catalog used for 
calibration. 

The 65 mJybeam^^ noise levels in the SPAM images 
match the lowest noise levels of the more than 500 im- 
ages that define the VLSS survey. A potential reduc- 
tion of the average noise level from 100 mJybeam^^to 
65 mJybeam"^ for the full VLSS survey would significantly 



increase the number of source detections from '^70,000 to 
about 120,000 (an increase of ~ 75%), but also it would 
greatly enhance virtually every science goal. For example, 
using the radio luminosity function for high-luminosity ra- 
dio galaxies from Jarvis et al. I|200ip . the estimated number 
of detectable HzRGs in the VLSS would increase by 65%, 
but also the maximum rcdshift would increase. For a lumi- 
nous radio galaxy with luminosity of 2 x 10^* WHz^^ sr^^ 
at 74 MHz, the redshift limit would rise from z = 5.7 to 
z = 6.8. Another example is the detection and study of 
cluster radio halos. Using available halo population models 
(Enfihn & Rottgcring lMI^ Cassano et al. lMIB)) . the antic- 
ipated noise reduction would roughly double the number of 
detectable halo systems. 

For the VLSS, the estimated theoretical thermal noise 
level of 35 mJybeam~^is still a factor of two lower than 
the average background noise level of ~ 65 mJybeam^^ in 
the SPAM images. From inspection of the SPAM images we 
cannot identify an obvious single cause for this. Therefore, 
similar to Cohen et al. 120071 we expect the remaining excess 
noise to be the combined result of several different causes. 
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including residual ionospheric phase errors after SPAM cal- 
ibration, but also residual RFI, collective sidelobe noise 
from many non-deconvolved sources (too faint or outside 
the FoV) and variable source amplitude errors (e.g., due 
to pointing errors and non-circular antenna beam patterns; 
see Bhatnagar et al. I2008P . 

The SPAM test results indicate that the ionospheric cal- 
ibration accuracy may be further improved. The typical 
model fit RMS phase error per antenna of ^ 20 — 30 de- 
grees for real data sets is much larger than the 3 degrees for 
the noiseless simulated data set. There are several possible 
sources of error, either in the peeling phase corrections or 
the ionospheric phase model. Noise in the visibilities (ei- 
ther thermal or non-thermal), contamination from other 
sources, inaccuracies in the peeling source model and un- 
dersampling of the fastest phase fluctuations are factors 
that degrade the accuracy of peeling. Also, the ionospheric 
phase screen model may be a poor representation of re- 
ality, either because it is incomplete (e.g., absence of ver- 
tical structure) or the fixed model parameters arc chosen 
poorly (e.g., screen height, spectral index of phase fluctu- 
ations). Several of these issues will be addressed in future 
work (Section [HI) . 

The potential problems with the peeling technique 
raises the question whether one should use alterna- 
tive methods. Apart from the precautions described in 
Section [3731 we have found little means to improve the ac- 
curacy of the peeling process for single sources any further. 
One unexplored option is to peel sources in groups, e.g. 
identify isoplanatic patches of sky with a large enough to- 
tal flux from multiple sources. Two possible alternatives 
approaches to peeling are: (i) simultaneous self-calibration 
towards multiple sources in the FoV, or (ii) fitting the iono- 
sphere model directly to the visibilities rather than using 
peeling as an intermediate step. Although these alternative 
approaches have not been tested by us in practise, we an- 
ticipate little improvement over our current accuracy. Van 
der Tol, Jeffs & van der Veen (|2007p show that, theoret- 
ically, iterative peeling converges to the same solution as 
simultaneous self-calibration. A direct fit of the ionosphere 
model to the visibilities is, similar to self-calibration, biased 
towards accurate solving in the direction of the apparently 
strongest source in the FoV. Although not conclusive for 
this approach, tests with SPAM show that using even a 
moderate flux-based weighting into the ionospheric phase 
model fitting against peeling phase corrections introduces a 
strong bias towards the brightest source, while calibration 
accuracy towards other peeled sources degrades severely. 

For the existing and future large low-frequency radio in- 
terferometer arrays like VLA-A, GMRT, LOFAR, LWA and 
SKA, the need for a direction-dependent ionospheric cali- 
bration method is evident. Based on the results presented 
in this paper, it is difficult to draw quantitative conclusions 
on the achievable calibration accuracy for these arrays. If 
a SPAM-like calibration algorithm is to be used in a very 
high signal-to-noise observing regime under quiet to mod- 
erate ionospheric conditions, it seems likely that residual 
RMS phase errors in the order of a few degrees could be 
achieved, comparable to the SPAM results on the simulated 
VLSS data set. 

When relying on the array itself to provide the necessary 
measurements to constrain ionospheric correction models, 
ionospheric calibration requires an array layout and sen- 
sitivity that allows for sampling the ionsphere over the 



array at the relevant spatial scales and time resolution. 
The spatial sampling is determined by the instantaneous 
pierce point distribution (or more general, the distribution 
of lines-of-sight through the ionosphere), which depends on 
the array layout and the detectable calibrator constella- 
tion. For future design of low-frequency arrays, it is recom- 
mended to optimize the array layout not just for scientific 
arguments (in general, centrally dense and sparse outside 
for good UV coverage), but also for ionospheric calibrability 
(in general, both uniform and randomized). 

6. Future Work 

To test the robustness and limitations of the method, it 
is necessary to apply SPAM calibration on a wide variety 
of data sets at different (low) frequencies, obtained with 
different arrays under different ionospheric conditions. Our 
highest priority is to test SPAM on observations from the 
largest existing LF arrays; the VLA in A-configuration and 
the GMRT. Data for these tests have been obtained and 
tests are currently in progress. One important possible lim- 
itation is the use of a 2-dimensional phase screen to repre- 
sent the ionosphere. We plan to expand the SPAM model by 
including multiple screens at different heights and compare 
the resulting image properties against the current single 
screen model. 

Another limitation of the current implementation is the 
absence of restrictions on the time behaviour of the model. 
Antenna-based peeling phases clearly show a coherent tem- 
poral behaviour, which is likely to exist for physical reasons. 
This could be used to reduce the number of required model 
parameters and suppress the noise propagation from the 
peeling solutions. We are currently investigating the pos- 
sibilities of forcing the SPAM model to be continuous in 
time. 

Several of the authors of this article are currently in- 
volved in setting up a simulation framework in which one 
has full control over the sky emission, ionospheric behaviour 
and array characteristics when generating artificial low fre- 
quency observations. Like in the test case on simulated data 
presented in Section |4l this allows for direct and quantita- 
tive comparison between the distorting ionosphere model 
and the recovered ionospheric phase model by SPAM. We 
plan to use this setup to further test optimize SPAM cali- 
bration for a broad range of ionospheric conditions. 
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Appendix A: Derivation and interpolation of the 
KL base vectors 

This Section contains an outline of the derivation and in- 
terpolation of the Karhunen-Loeve (KL) base vectors that 



are used to describe the ionospheric phase screen in SPAM. 
The KL technique is adopted from the work by van der 
Tol & van dcr Veen ( |20(J7i) • For a given stochastic model of 
spatial electron density fluctuations in the ionosphere, the 
differential phase rotation on rays of passing radio waves 
can be described by a (zero mean) isotropic phase screen 
(j){p) with a given spatial covariance function 



(A.l) 



where (...) denotes the expected value and r = \r\ is the 
length of vector r. In Kolmogorov turbulence theory, a 
phase structure function is defined as (Equation [5]) 



(A.2) 



The structure function and the covariance function are re- 
lated through 
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(A.3) 



where C^^{0) = a'^ is the phase variance. The phase struc- 
ture function of a plane wave that passed through a tur- 
bulent layer is found to behave as a power-law over a large 
range of spatial scales (see Section [^?T|) 



(A.4) 



where Tq is a measure of the scale size of phase fluctuations 
and 7 is the power-law slope. 

Because the domain of p is a limited set of P discrete 
ionspheric pierce points p, we switch to matrix notation. 
The elements of the [P x P) phase covariance matrix are 
given by 



(A.5) 



with Pi^Pj e {p} and = \p^ ~ Pj\- The symmetric co- 
variance matrix can be decomposed into 



(A.6) 



where the columns of (P x P) matrix U contain the or- 
thonormal eigenvectors of C^^. U"'" is the transpose of U 
and the {P x P) diagonal matrix A contains the eigenval- 
ues. The eigenvectors in the columns of U are a suitable set 
of base vectors to describe the phase screen at the pierce 
points {p}. The eigenvalues on the diagonal of A are a mea- 
sure of the variance of the eigenvector coefficients. Reducing 
the fitting order is an arbitrary, but necessary step. An op- 
timal subset of eigenvectors is determined by selecting only 
those with the largest eigenvalues, as the coefficients with 
the highest variances are the most significant in the mod- 
eling problem. 

When the fitting order is reduced from P to P < P, 
we are left with a subset of eigenvectors in the columns 
of {P x P) matrix U, and eigenvalues on the diagonal of 
(P X P) matrix A, for which 

C^^^UAU^. (A.7) 

The phase screen at the pierce point locations can be ap- 
proximated by 



(A., 
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where we have denoted the phases and eigenvcetor coef- 
ficients in matrix notation as $ (P x 1) and q (P x 1), 
respectively, q is unknown and needs to be solved for by 
use of a non-linear least squares method using Equation [211 
When q is determined, the phase screen can be evaluated 
at arbitrary pierce point locations {p} through Kriging in- 
terpolation fMatheron ll973p : 



* = C - C-]Uq 



(A.9) 



where C - is the (P x P) covariance matrix between {»} 

ipip 

and {p}. Inversion of the full can be approximated 
using 

C7i 



(A.IO) 



The elements of the covariance matrix can be calculated 
using Eauations lA.3l and lA.4l For our application, the abso- 
lute value of 7'g is not relevant because we only require the 
relative eigenvalues for the order reduction. The elements 
of the (P X P) structure matrix are given by 



D 



00[*>JJ = (A.ll) 

The relation between the structure matrix and the covari- 
ance matrix is given by 

1. 



(A.12) 



with 1 a (P X 1) vector containing ones. The phase variance 
terms cr^ are removed from the equations by explicitly re- 
moving the mean phase from the individual phases through 
the substitution 



(I- -^11^ 



(A.13) 



where I is the (P x P) identity matrix. Applying this sub- 
stitution to the covariance matrix yields 



1 H 
-11^ 

P 

-11^ 

P 



1 n 
— 11^ 

P 



- -11^ 

P 



(A.14) 



The CT? terms have dropped because of the properties 



(I-lllT)l = 0, 



^(I-lllT)=0 



(A.15) 



For Kriging interpolation, a similar substitution is per- 
formed as in Equation lA. 131 



(A.16) 



with i a (P X 1) vector containing ones. The covariance 
matrix between {p} and {p} (Equation IA.9P now becomes 



2 </"/> 



1 ~ n 
-11^ 

P 



1 H 

p'' 



(A.17) 



where D^^ is the (P x P) structure matrix between {p} 
and {p}, calculated using Equations I A. 41 and I A. 1 II 



